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Simulations of a supersonic recessed-cavity flow are performed using a hybrid large- 
eddy/Reynolds-averaged simulation approach utilizing an inflow turbulence recycling 
procedure and hybridized inviscid flux scheme. Calorically perfect air enters the three- 
dimensional domain at a free stream Mach number of 2.92. Simulations are performed to 
assess grid sensitivity of the solution, efficacy of the turbulence recycling, and effect of the 
shock sensor used with the hybridized inviscid flux scheme. Analysis of the turbulent 
boundary layer upstream of the rearward-facing step for each case indicates excellent 
agreement with theoretical predictions. Mean velocity and pressure results are compared to 
Reynolds-averaged simulations and experimental data for each case, and these comparisons 
indicate good agreement on the finest grid. Simulations are repeated on a coarsened grid, 
and results indicate strong grid density sensitivity. The effect of turbulence recycling on the 
solution is illustrated by performing coarse grid simulations with and without inflow 
turbulence recycling. Two shock sensors, one of Ducros and one of Larsson, are assessed for 
use with the hybridized inviscid flux reconstruction scheme. 
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Nomenclature 

Klebanoff function constant 
turbulence model constant 
cavity depth 

distance to nearest solid surface 
cell interface flux 
von Karman constant 
turbulence kinetic energy 
shock sensor constant i 
pressure 

primitive variable vector 
time 

temperature 

streamwise velocity component 
transverse velocity component 
Cartesian coordinates 

constant controlling the LES/RAS blending 
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Kronecker delta function 
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boundary layer thickness 
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density 
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LES/RAS blending function 

r ss 
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shock sensor coefficient 
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kinematic molecular viscosity 

CO 
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turbulence frequency 
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Taylor micro-scale 

inf 
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free- stream value 
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= 

spatial coordinates 

Kleb 
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Klebanoff 

res 

= 

resolved 

SGS 
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subgrid-scale 

HY 

= 

hybrid 

SR 

= 

symmetric reconstruction 

UP 

= 

upwind 


= 

fluctuation 

* 

= 

normal to a surface 

+ 

= 

wall unit 


I. Introduction 

In order to address some of the fundamental limitations of Reynolds-averaged simulation (RAS) approaches while also 
exploiting the strengths of large-eddy simulation (LES) approaches, researchers are more frequently using a blended 
LES/RAS approach. This is especially true for use with high Reynolds number supersonic wall-bounded flows, 
which exhibit inherent unsteadiness, turbulent shock-boundary layer interactions, and other complex flow 
phenomena. Whereas the RAS equations require modeling of turbulence at all length scales, the LES equations only 
require modeling of turbulent motions of a length scale smaller than a characteristic filter width; the largest turbulent 
motions are resolved naturally through the solution of the filtered governing equations. Since these larger turbulent 
structures are responsible for the majority of mass, momentum, and energy transport, LES simulations are well- 
suited for high Reynolds number flows. Unfortunately, the grid densities required for pure LES simulations of wall- 
bounded flows often render the approach computationally infeasible for problems of practical importance. 

In this study, a hybrid LES/RAS approach is used in conjunction with a turbulence recycling procedure to 
simulate the Settles, et. al. 1,2,3 supersonic flow over a rearward facing step and subsequent reattachment along an 
inclined wall. The recycling process allows for the introduction of resolved turbulent content in an equilibrium 
turbulent boundary layer without having to simulate the entire development of the boundary layer. The chosen 
simulation framework is assessed in its ability to predict mean flow properties within the approach turbulent 
boundary layer, through the detached shear layer, and along the inclined wall. Assessments of the LES/RAS 
simulations are made using steady-state RAS results and experimental data. 

This study serves a dual-purpose; in addition to validating the application and extension of the hybrid 
LES/RAS methodology, these simulations assess several recent developments to the Viscous Upwind aLgorithm for 
Complex flow ANalysis (VULCAN) computational fluid dynamics (CFD) code. 4 These recent developments 
include reduced-dissipation numerics, a hybridized flux reconstruction scheme, and inflow turbulence recycling. 
Specific objectives of the current study include: assessment of the applicability of the hybrid LES/RAS method for 
cavity flows, isolation of the effect of inflow turbulence recycling on the flowfield statistics, determination of any 
grid sensitivity of the solution, and evaluation of shock sensor choice for use with the hybridized inviscid flux 
reconstruction scheme. These objectives are achieved through the study of six distinct cases that are parameterized 
by solution method, inflow treatment, grid size, and shock sensor and are summarized in Table 1. Elsewhere in this 
paper, each of the simulation cases will be referenced using the case ID stated in Table 1. For example, references 
will be made using “case R35D”, rather than “case one”, where R, 35, and D indicate recycling, 35 million grid 
cells, and Ducros sensor, respectively. 
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Table 1. Simulation cases parameterized by inflow treatment, grid size, and shock sensor. 


Case 

Case ID 

Solution Method 

Inflow Treatment 

Grid Cells 

Shock Sensor 

1 

R35D 

FES/RAS 

Recycling 

35e6 

Ducros 

2 

RR8D 

FES/RAS 

Recycling 

8.8e6 

Ducros 

3 

NR8D 

FES/RAS 

No recycling 

8.8e6 

Ducros 

4 

R35F 

FES/RAS 

Recycling 

35e6 

Ears son 

5 

R8F 

FES/RAS 

Recycling 

8.8e6 

Ears son 

6 

NR8RAS 

RAS 

No recycling 

8.8e6 

None 


In the following section, the original wind-tunnel experiment and major flow features are described. In 
section three, the hybrid LES/RAS formulation and associated turbulence closure models are described. Details of 
the numerical formulation are included in section four. After describing the grids and computational executions in 
section five, the results and accompanying discussion are presented in section six. Conclusions and summarizing 
comments are included in the final section. 


II. Physical Flow 

The flow chosen for the current investigation was first studied by Settles, et. al. in the early 1980s. 1 ’ 2,3 The 
original experiment was performed in the Princeton University 20cm by 20cm High Reynolds Number Supersonic 
Wind Tunnel and consisted of Mach 2.92 air entering the tunnel test section at a stagnation temperature of 258K, a 
stagnation pressure of 0.69MPa, and a free stream unit Reynolds number of 6.7x1 0 7 /m. The test article, presented 
in Fig. 1, was designed such that a turbulent boundary layer developed along a flat plate before separating over a 
sharp corner. The free shear layer formed over the subsequent cavity and reattached along a wall of length 18cm 
inclined at 20°. The position of the inclined wall was adjusted before taking measurements until no change in 
pressure or flow direction occurred due to the separating boundary layer. This procedure was necessary due to the 
high sensitivity of the shear layer location to the pressure through the test region. The cavity depth was 2.54cm, and 
the model was unaffected by the tunnel sidewall boundary layers due to its positioning away from the walls. In 
order to further promote two-dimensionality, aerodynamic fences lined the edges of the inclined wall. Mean flow 
measurements were made using a traversing set of pitot and static pressure probes, and these measurements, in 
conjunction with an assumption of constant total temperature, enabled the researchers to calculate Mach number, 
velocity, density, and other properties. A spark shadowgraph enabled observation of major flow features. 

The original experiment was designed to illustrate the complex fluid phenomena comprising the reattachment of 
a free shear layer in a supersonic flow. The effect of turbulence on the flowfield is prominent, as documented in the 
original review of the experiment and in several subsequent numerical investigations. 2,5 Major flow features 
observed in the experiment are illustrated in Fig. 2. After the approach boundary layer detaches at the step corner, a 
free shear layer forms over the cavity. As the shear layer reattaches at the inclined wall, fluid is entrained into the 
cavity, thereby forming a recirculation zone. The inherent unsteadiness of the shear layer alters the location of 
reattachment and subsequently affects the fluid entrainment. Together, the shear layer unsteadiness, altered 
attachment point, and recirculation zone form a feedback loop that gives rise to a low-frequency flapping motion of 
the shear layer. Furthermore, as the shear layer compresses and reattaches, an oblique compression shock is formed 
off the inclined wall. 



Figure 1. Experimental setup with one aerodynamic fence removed and the computational domain footprint 
signified by dashed lines. 
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M=2.92 



Figure 2. Major flow features observed during the experiment. 


III. Hybrid LES/RAS Formulation 

In the current hybrid LES/RAS approach, blending of the solution between LES and RAS throughout the 
physical domain is achieved via blending of the closure models. The current implementation in VULCAN is 
designed such that the choice of RAS and SGS closure models is delegated to the user, thereby allowing significant 
flexibility. An overview of the detailed description of this formulation found in Ref. 6 is presented here. 

VULCAN solves the Navier- Stokes equations using a cell-centered finite- volume scheme and an upwind 
formulation for calorically or thermally perfect gases. Various one-equation and two-equation turbulence models 
are available for treating the turbulence closure. The code has been used extensively for simulating high-speed 
reacting flows in its RAS capacity with much success. However, by its nature, the RAS formulation is inadequate 
for the study of strongly unsteady flows, and the turbulence models used to close the RAS equations are often 
questionable for complex flowfields of engineering interest. The current hybrid LES/RAS capability enables 
VULCAN to be used more reliably in highly unsteady applications, in which large-scale turbulence and transients 
play an important role in simulating and understanding the evolution of a flowfield. 

The blending of the LES and RAS approaches is accomplished using a flow- and grid-dependent blending of the 
chosen turbulence closure models. In this work, the Menter-BSL turbulence model 7 is used to close the RAS 
equations. This model requires the solution of two additional equations — a transport equation for turbulent kinetic 
energy and a transport equation for turbulence frequency. To model the subgrid-scale (SGS) turbulence for use in 
LES regions, the one-equation model of Yoshizawa and Horiuti 8 is used. This model requires the solution of a 
transport equation for the SGS turbulence kinetic energy (TKE). A detailed description of the turbulence models 
and the complete set of hybridized governing equations are found in Ref. 5. 

The LES/RAS blending function is designed such that the RAS formulation transitions to the LES formulation in 
the outer region of the boundary layer, near the upper edge of the log layer. As previous studies have indicated 9 , 
targeting the transition for this region at the edge of the log layer within the turbulent boundary layer yields the most 
resolved turbulence without adversely affecting the universality of the boundary layer. The blending from LES to 
RAS is done smoothly using the following function 10 : 


r = i(l — tanft [5(^„ 2 -l) -</>]), „=^ 


( 1 ) 


where the Taylor microscale is defined as: 


X = 



( 2 ) 


The blending function input parameter, a , is set by the user such that r = 0.99 at the targeted average LES/RAS 
transition point, and is set to atanh(0.98) to fix the balancing position, where k r] 2 = to T = 0.99. Performing 
the procedure outlined in Ref. 1 1 for the chosen flow yields a value of a = 32.8 for the turbulent approach boundary 
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layer, which forces the switch from RAS to LES at the outer edge of the log layer (corresponding to a d + of about 

1080 ). 


IV. Numerical Formulation 

The simulations performed for this investigation are an early application of several recent numerical 
developments to the VULCAN CFD code. VULCAN was originally designed to serve as a robust, upwind, RAS 
solver for application to high-speed reacting and non-reacting flows. In recent years, VULCAN has been expanded 
to include LES and hybrid LES/RAS capabilities, which drove an investment 
in various higher-order numerics. Several of these methods are used in the current study and are described in detail 
below. 


A. Low-Dissipation Numerics and Hybridized Inviscid Flux Scheme 

In order to reduce numerical dissipation in smooth flow regions, while maintaining the dissipation necessary to 
capture discontinuities elsewhere, a hybridized flux scheme is used to reconstruct interface inviscid fluxes. This 
hybridization is implemented as the linear combination of a low-dissipation, symmetrically-reconstructed inviscid 
flux and a conventional upwind-biased formulation. The low-dissipation flux formulation is described by Gieseking 
et al ., 12 in which the left and right states are reconstructed to higher-order accuracy at the cell interfaces using a four- 
point symmetric interpolation on the primitive variables of the form: 

$“+l = C+i = ii Mi + 3<+i) - 7i Mi+2 - 9i-i) < 3 > 

The low-dissipation interface flux is computed using a conventional average of the left and right fluxes according to: 



(4) 


This low-dissipation flux, while accurate in regions of smooth flow, is incapable of resolving shock discontinuities 
and therefore must be combined with a more conventional dissipative flux. 

To construct the necessary dissipative flux, the piecewise-parabolic method of Collela and Woodward 13 is first 
used to reconstruct the primitive state variables at the left and right sides of the cell interface. A minmod-type 
limiting procedure 9 is used as a means of capturing discontinuities. With the higher-order states computed at the 
interfaces, the Harten, Lax, van Leer, and Contact (HLLC) scheme is then used to compute the dissipative flux . 14,15 
The HLLC scheme is a successor to the Harten, Lax, and van Leer (HLL) scheme and is designed to include 
treatment of the missing contact surface. 

With the low-dissipation and dissipative fluxes computed, the final hybridized inviscid interface flux is 
computed according to: 


f h \ = r ss F ut [ + (i - f ss )f sr 1 

i+— ^ i+- i+- 

2 2 2 


(5) 


where L S s is a flow-dependent scaling factor and is referred to as a “shock sensor” due to its functional purpose. As 
briefly discussed before, application of LES to high-Mach number flows requires dissipative schemes for capturing 
shocks and low-dissipation schemes for resolving vortical motions. A shock sensing switch is thus used to 
dynamically adjust the solver for conditions near to and away from a shock. As illustrated above, the shock sensor 
is used to scale the respective flux components comprising the hybridized inviscid flux reconstruction. In this study, 
two methods for sensing shocks are used: the Ducros sensor 16 and the approach of Larsson . 17 The Ducros sensor is 
computed as: 


r _(V-U) 2 / 

lss ~ / (V • it ) 2 + (V x u ) 2 + 


( 6 ) 


where is a small number preventing a divide-by-zero operation when flow gradients are not present. The Larsson 
sensor is calculated according to: 
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(7) 


l ss 


= 1 ? 


0 for 9 < 1 
for 0 > 1 


9 = 


-(V-u) 


max (l sw |Vxu|,L ss>2 £) 


( 8 ) 


where a is the sound speed and ds = \/cell volume. The constants L ss l and L ss2 were set to 5.0 and 0.05, 
respectively, for the current study. Furthermore, the lower bound for each of these sensors was modified to provide 
a small level of background dissipation. Rather than allowing it to range from 0.0 to 1.0, each shock sensor is 
bounded by 0.2 and 1.0, thereby preventing numerical instabilities arising within the shear layer reattachment 
region. 

B. Turbulence Recycling Procedure 

In order to introduce coherent turbulent structures at the inflow, turbulent fluctuations and relevant 
turbulence properties are recycled from a plane one cavity-depth upstream of the step. Turbulent fluctuations are 
extracted from the donor plane by subtracting time- and span-averaged mean flow data from the instantaneous flow 
data. These fluctuations are scaled according to boundary-layer similarity laws before being superimposed to a 
mean flow RAS profile. The general procedure used for this process is described by Edwards and others, 9 ’ 1819 and 
several important characteristics of the procedure are described here. 

In order to prevent excessive turbulent energy accumulation in the outer regions of the boundary layer, an 
intermittency factor is applied to the recycled fluctuations. This intermittency factor is applied to density, 
temperature, and velocity components according to Eq. (9) and to turbulent kinetic energy, turbulence frequency, 
and the hybrid switching function according to Eq. (10). 


Qinflow RrAS ^KlebRrecyc R P)T,Ui 

Qinflow — (1 — ^Kleb^RAS ^ KlebQrecyc 9 — k,(i),r 

,61 


Fltleb — 


1 + 


( d /c Kleb s) _ 


-1 


(9) 

( 10 ) 

( 11 ) 


In the preceding equations, q is a placeholder for the variables referenced above. The intermittency factor, F k i eb , is a 
Klebanoff-type intermittency function, 20 where C k i e b = 1-10. This modification ensures that the RAS free stream 
inflow properties are not modified outside the boundary layer by the recycling procedure. 

A further modification of the general recycling procedure is to consider pressure fluctuations as they relate 
to the calculation of density fluctuations. Following the Morkovin hypothesis, 20 which assumes that total 
temperature fluctuations are negligible, the temperature fluctuations are related to the velocity fluctuations. The 
recycled temperature fluctuations are then limited to prevent their magnitude from exceeding that predicted by 
Morkovin’ s hypothesis. Using these limited temperature fluctuations, a provisional pressure fluctuation is 
calculated, which is also limited according to some user-defined multiple of the mean static pressure (chosen as 0.02 
in the present simulations). The limited recycled temperature and pressure fluctuations are used in conjunction with 
the mean temperature and pressure to compute corrected density fluctuations. 


V. Computational Execution 

Simulations are performed on a structured three-dimensional grid consisting of approximately 35 million cells 
for the fine level and approximately 8.8 million cells for the coarsened level. The coarse grid was constructed by 
coarsening the fine grid once in each of the streamwise and spanwise directions, where coarsening is defined as 
removing every other grid point in a given direction. The grid is composed of four primary blocks, as shown in Fig. 
3. The two upper blocks are connected to the lower blocks via non-CO continuous patches. Similarly, the lower left 
block is attached to the lower right block by a non-CO continuous patch where only the spanwise coordinates are 
mismatched (the downstream spacing is twice that of the upstream spacing). For the fine grid, the I-J-K 
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(streamwise-transverse-spanwise or x-y-z) dimensions of the blocks adjacent to the wall directly upstream and 
downstream of the step are 525x121x301 and 453x241x151, respectively; dimensions of the upper left and right 
blocks are 29x45x15 and 153x45x15, respectively. The patches connecting the lower and upper blocks are 
positioned far enough from the domain of interest to ensure that resolved boundary layer and free shear layer eddies 
are not convected into them. The physical dimensions in the x, y, and z directions of the computational domain are 
28.6cm, 22.5cm, and 3.81cm, respectively. The flat plate extended four step heights upstream of the step, which is a 
factor of two further upstream than in previous studies 5,21 to allow for recycling of turbulent fluctuations, and the 
grid was clustered to provide a nominal y + of 1.0 adjacent to the surface. The streamwise and spanwise grid spacing 
was set to l/20 th of the boundary layer thickness. Downstream of the step, the grid clustering required by the 
approach flow boundary layer was relaxed toward a more isotropic grid for the resulting shear layer. The maximum 
x, y, and z spacing in this region was 2.5%, 5.0%, and 1.0%, respectively, of the cavity depth (2.54cm). The cell 
aspect ratio for the LES regions was no greater than approximately 20.0. The three-dimensional grid was split into 
1980 sub-blocks to allow runtime parallelization using Message Passing Interface (MPI) on NASA Advanced 
Supercomputing resources. Each simulation performed on the fine grid required approximately 240 hours of wall 
clock time on 512 processors. The same simulations performed on the coarse grid required about 60 hours of wall 
clock time on 512 processors. 

Inflow variable profiles and flowfield initializations for the LES/RAS cases were obtained using the converged 
RAS results of case NR8RAS. For this case, the Menter k-co turbulence model 7 was used with a compressibility 
correction 22 , and the free stream flow was allowed to develop until the momentum thickness one cavity-depth 
upstream of the step matched experimental measurements. A significant spread in solutions was obtained using 
various turbulence models for the preliminary RAS simulation, and the RAS results used in this study represent 
those in best agreement with experimental data. Nominal inflow conditions included: static pressure of 21.2kPa, 
Mach number of 2.92, and temperature of 92.4K. For all LES/RAS cases, homogeneous turbulence was assumed in 
the spanwise direction allowing the use of periodic boundary conditions at both boundaries. For the fine grid cases, 
all walls were assumed adiabatic and no slip. For the coarse grid cases, all walls except for the vertical cavity wall 
were assumed adiabatic and no slip; wall functions were applied at the vertical cavity wall due to excessively high 
y+ values observed there during runtime. For all runs, an extrapolation condition was applied at the exit plane, and 
a characteristic condition was enforced at the far-field plane. The inflow profiles obtained using case NR8RAS 
were applied at the inflow plane for each subsequent case. For cases in which turbulence recycling was enabled, 
turbulent fluctuations were superimposed onto the initial RAS solution upstream of the cavity step. These turbulent 
fluctuations were extracted from a previously-simulated fully-resolved (via LES) turbulent boundary layer on a flat 
plate. This initialization procedure provides the turbulence necessary for effective recycling without necessitating 
the natural development of such structures. 



x 


Figure 3. Isometric view of the three-dimensional domain coarsened twice in the streamwise and transverse 
directions and once in the spanwise dimension for visual clarity and to emphasis grid design. 
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The governing equations were integrated using a dual-time-stepping approach, 23 in which the diagonalized 
approximate factorization (DAF) method 24 is used for integration in pseudo time with a maximum Courant- 
Friedrichs-Lewy (CFL) number of 100.0, and a second-order three-point backwards finite difference approximation 
was used for integration in physical time. Convergence was defined as achieving a residual reduction of 2.5 orders- 
of-magnitude or as reaching a limiting number of 15 sub-iterations. The physical time step was set to 0.1 ps. 
Following initialization, the flowfield was evolved for a minimum of 10 characteristic flow-through times before 
taking statistics. A characteristic flow-through time was defined as the time it takes a particle in the free stream to 
travel from the corner of the step to the point of reattachment of the free shear layer along the inclined wall — 
approximately 0.2ms. Once a statistically stationary state was achieved, flowfield properties were ensemble- 
averaged in time at a rate of 0.5ps for periods of 5, 10, and 15 flow-through times and subsequently averaged 
spatially in the homogeneous spanwise direction. The RAS turbulence model constants were set according to Ref. 
7, and SGS turbulence model constants were set as suggested in Ref. 8. 

VI. Results and Discussion 

Fifteen characteristic flow-through times were necessary to achieve reliable statistics. As illustrated in Fig. 4, 
the simulations reproduced the major flow phenomena observed in the experiment; the supersonic flow upstream of 
the cavity detaches at the step corner, thereby creating a free shear layer over the cavity. This free shear layer 
attaches along the inclined wall and interacts with an oblique shock front standing off the inclined wall. Fluid is 
entrained into the cavity near the location of shear layer attachment, and this fluid creates a recirculation zone which 
drives an unsteady low-frequency motion of the free shear layer. Note that in Fig. 4, annotations are included that 
indicate major flow features and identify each image’s relative location in time, which is given as the value of 
elapsed time relative to an arbitrary origin, At, divided by the characteristic flow-through time, Ax. In the following 
sections, results are presented and organized according to study objective, including: assessment of solution 
accuracy via comparison to experiment, assessment of grid sensitivity, isolation of inflow turbulence recycling 
effect, and determination of shock sensor effect on solution accuracy. 

A. Comparison with Experiment 

Before comparing results to experimental data, the boundary layer was confirmed to have developed properly 
upstream of the step. This confirmation was necessary to ensure that the constant, a, required for the recycling 
procedure was chosen appropriately and that the recycling procedure itself did not adversely affect the inflow 
boundary layer. For both the fine grid and coarse grid cases, boundary layer velocity profiles at several planes 
upstream of the step were compared to theoretical predictions based on Cole’s Law-of-the-Wall formulation and the 
original RAS solutions under the van Driest transformation described in Ref. 9. These three planes included the 
inflow plane, the donor recycling plane, and a plane at the step corner. For each case and at each plane, the 
averaged LES/RAS results agreed well with both the theoretical log-layer predictions and with RAS data. The 
comparisons for cases R35D and R8D at the recycling plane are included in Fig. 5 to demonstrate the level of 
agreement. In each case and at each plane, the LES/RAS simulations reproduced the slope and location of the log- 
layer well, and the blending function anchored the mean LES/RAS transition point near the upper edge of the log- 
layer. Thus, the LES/RAS approach and recycling procedure provide proper development of the boundary layer 
upstream of the step. 
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Figure 4. Instantaneous snapshots of velocity magnitude at the spanwise two-dimensional centerline plane for 
case R35D with recycle plane, major flow features, and data extraction locations marked. 
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Figure 5. Comparisons of velocity and blending function through the approach boundary layer for cases 
R35D and R8D at the recycling plane. 

With the upstream boundary layer validated against theoretical predictions, downstream data obtained from the 
simulations are compared with that of experiments. Case R35D — the fine grid case using the Ducros sensor — is 
designated the baseline case. In Figs. 6a through 6d, velocity and pressure results are compared to available 
experimental data and to the compressibility-corrected RAS data for case NR8RAS. Note that a variable having a 
superscripted * indicates a value normal to a solid surface, and all spatial dimensions are normalized by the cavity 
depth, D. Overall, the LES/RAS results agree quite well with the experimental data after shear layer reattachment, 
and a mismatch is noted in the shear layer profiles. The deviation apparent in these results is an under-prediction of 
the vertical location of the shear layer, which is most-clearly visible in the shear layer velocity profiles. This 
deficiency also appears in the velocity and pressure profiles, normalized by free stream values, normal to the 
inclined wall near the point of reattachment (x /D=2. 7). The assumption of turbulent homogeneity in the spanwise- 
direction may be the cause for this issue, since previous efforts have shown that three-dimensional effects play a 
significant role in flow separation even when extensive steps are taken in the experiments to ensure some level of 
two-dimensionality. 25 In this experiment, flow fences were placed on either side of the compression ramp. 
However, the interaction of the reattaching shear layer with theses fences was not captured in the current 
simulations. It is also plausible that the width of the homogeneous direction was not large enough to adequately 
capture the largest scales associated with the recirculation region within the cavity. A future effort is planned to 
assess this possibility. Predictions of the normalized static pressure along the inclined wall are considerably better 
for LES/RAS, as well as predictions of the normalized velocity and pressure profiles at planes extending normal to 
the inclined wall. This result was expected given the fact that RAS models tend to severely over-predict the 
boundary layer recovery rate when free shear flows reattach to solid surfaces. 
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Figure 6. Comparisons of normalized velocity and pressure for cases R35D (blue dashed line) and NR8RAS 
(solid red line) to experiment (open triangles). 


Though the vertical location of the shear layer is under-predicted, the shear layer spreading rate is captured quite 
accurately by the LES/RAS data. This can be seen by translating the velocity curves through the shear layer 
vertically until the uppermost point of inflection of the simulation data intercepts the experimental data. Figure 7 
shows these vertically-translated curves, which demonstrate the difference in spreading rate predictions between 
RAS and LES/RAS. Assuming the shear layer forms a line across the cavity, the simulated shear layer for the 
baseline case is deflected about 1.8° downward from the experimental data. The good spreading rate agreement, 
when considered in light of the mismatch in vertical displacement of the shear layer, strengthens the previous 
discussions on the design differences between the experiment and simulations. Not only is it likely that the 
experiment’s aerodynamic fences only partially alleviate three-dimensional effects, but the limited physical extent of 
the computational domain’s spanwise direction may cause further deviation from experiment due to artificial 
limiting of the size of the largest spanwise flow features that can reside within the recirculation region. 
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Figure 7. Normalized velocity profiles through the shear layer translated vertically to demonstrate spreading 
rate agreement for cases R35D (blue dashed line) and NR8RAS (red solid line), as compared to experimental 
data (open triangles). 


B. Grid Sensitivity 

An assessment of the grid sensitivity was made by comparing results of cases R35D and R8D, as the only 
difference between these runs was the number of grid cells used to represent the physical domain. The grid used for 
case R8D was identical to that of case R35D, albeit coarsened once in the streamwise and spanwise directions. The 
transverse direction was not coarsened because the resulting grid would have given y+ values too large for integrate- 
to-the-wall simulations. Each case was evolved in time for ten flow-through times before statistics were collected 
for 5, 10, and 15 flow-through times. The grid sensitivity of the solution was found to be significant. The coarsened 
grid is unable to resolve adequate levels of turbulence, as compared to the fine grid. The level of resolved 
turbulence captured by each grid is significantly different, as illustrated in Fig. 8, where instantaneous Mach number 
contours are presented at a two-dimensional slice at the spanwise centerline. Note that while the turbulent 
fluctuations in the approach flow boundary layer are well-resolved for both cases, these fluctuations and 
corresponding turbulent motions are resolved to a lesser degree moving downstream of the step when the coarse grid 
is utilized. This observation is further illustrated through examination of the contributions of resolved TKE and 
modeled TKE downstream of the step and is demonstrated in Fig. 9, where line plots of the TKE contributions are 
presented at several planes. At each of these planes, a measure of the resolved turbulent content is approximated 
according to: 


%TKE res = 


TKE res . „ „ 

— x 100 

TKEsgS +TKE res 


( 12 ) 


Calculations using Eq. (12) are performed using peak-to-peak values from Fig. 9. Whereas in case R35D, the 
turbulence is 71%, 80%, and 85% resolved at each plane, respectively, the same is not true for case R8D. For case 
R8D, the turbulence is only 30%, 36%, and 50% resolved at each of the three planes, respectively. Typically for 
LES of high Reynolds number flows, the percentage of resolved turbulence is in the 80-90% range. 26,27,28 The 
dominance of modeled TKE leads to more reliance on the SGS model and, in this case, a larger disagreement with 
the available experimental data. This disagreement is highlighted in Fig. 10, where the results of case R8D are 
plotted with the results of case R35D and the experimental data through the shear layer. Note that the data for each 
case have been shifted, as described earlier, to highlight the difference in spreading rates. 
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Figure 8. Instantaneous Mach number snapshots for cases R35D (top) and R8D (bottom) at the spanwise 
centerline plane. 


x/D = 0.7 x/D = 1.1 x/D = 1.9 





Figure 9. TKE contributions at three planes downstream of the step for cases R35D and R8D. The blue 
dashed line and red solid line represent the modeled TKE for cases R35D and R8D, respectively, and the gray 
dash-dot-dot line and green gash-dot line represent the resolved TKE for cases R35D and R8D, respectively. 
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Figure 10. Normalized velocity results through the shear layer translated vertically for cases R35D (blue 
dashed line) and R8D (black solid line), as compared to experiment (open triangles). 


C. Effect of Inflow Turbulence Recycling 

The effect of inflow turbulence recycling on the flowfield was isolated by comparing the results of cases R8D 
and NR8D. For these two cases, the only differentiating factor is the inflow turbulence recycling, which is enabled 
in the first and disabled in the second case, respectively. Specifically, for case NR8D, the RAS solution from case 
NR8RAS was used to specify the inflow conditions at the comer of the step; the LES/RAS governing equations 
were then solved in space and time for the flowfield downstream of the step. This treatment is similar to how a 
detached eddy simulation approach would proceed, where attached boundary layer regions are treated as pure RAS 
and only detached regions are treated with LES. Accounting for resolved fluctuations in the approach boundary 
layer had little effect on the shear layer development as demonstrated by the results of cases R8D and NR8D in Fig. 
11. The shifted velocity profiles through the shear layer are in relatively good agreement between the two cases. 
Keeping in mind that the grid used in cases R8D and NR8D inadequately resolved turbulence downstream of the 
step, the current results indicate the inflow turbulence recycling had little effect on the shear layer development for 
the current flow. 


x/D=-1 .0 


x/D=1.5 x/D=2.5 x/D=3.5 




Figure 11. Normalized velocity results through the shear layer for cases R8D (black solid line) and NR8D 


(magenta short-dash line), as compared to experiment (open triangles). 
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D. Effect of Shock Sensor 

In order to determine the effect that choice of shock sensor may have on the solution, cases R35D and R8D were 
repeated using the Larsson shock sensor in place of the Ducros sensor. The Ducros sensor was considered the 
baseline shock sensor for purposes of this work. In their implementation, these shock sensors compare the 
magnitude of the divergence and curl of velocity; however, the Larsson sensor uses a slightly different formulation, 
which includes logic that toggles the interface flux approximation to either fully-upwind or fully- symmetric. 
Conversely, the Ducros switch blends the two independent interface flux approximations in proportion to the 
magnitude of the normalized sensor coefficient. 

Instantaneous snapshots of the shock sensor for cases R8D and R8L are presented in Fig. 12. Note that the 
magnitude of the shock sensor differs significantly across the flowfield between the two cases. As explained 
previously, the Ducros sensor takes on values between the upper and lower bounds of the switch throughout most of 
the domain, unlike the Larsson sensor, which forms more distinct regions limited to the upper and lower sensor 
bounds. Though the sensor distributions are substantially different, the results of cases R35D and R35L, and of 
cases R8D and R8L, are in relatively good agreement. Normalized velocity and pressure profiles, through the shear 
layer and along the inclined wall, are presented in Fig. 13 for cases R35D and R35L. The velocity profiles in Fig. 
15a have been shifted to compare shear layer spreading rate and illustrate good agreement. Assuming a line 
connects the curve at x/D=3.5 to the comer of the step, the angular deviation between the two shear layers is about 
0.85°. The pressure profiles plotted in Fig. 15b exhibit a small difference, likely attributed to the vertical 
displacement of the shear layer. In case R35L, the shear layer intersects the vertical wall slightly lower than in case 
R35D, and this causes a strengthened oblique shock off the inclined wall. This drives a higher rise in pressure along 
the inclined wall and alters the velocity and pressure profiles normal to the inclined wall at the first data station, 
x*/D=2.7. Comparisons of cases R8D and R8L indicate even better agreement than those included in Fig. 15 and 
are therefore omitted here. Based on these comparisons, it is reasonable to conclude that the choice of shock sensor 
has little effect on the statistics for this flowfield. 




-2 0 2 x/D 4 6 8 


Figure 12. Instantaneous shock sensor contours for cases R8D (top) and R8L (bottom) at the spanwise 
centerline plane. 


15 

American Institute of Aeronautics and Astronautics 


x/D=-1.0 


x/D=1.5 


x/D=2.5 


x/D=3.5 


0.2 r 


0 . 1 - 


! 0 . 0 - 


- 0 . 1 - 


- 0.2 


a) 


a 

/A 


/a 

A 


1 0 


_L L_ 
1 0 

u/U 



x*/D=2.7 


x*/D=3.8 x*/D=4.6 x*/D=5.4 


x*/D=2.7 


x*/D=3.8 x*/D=4.6 x*/D=5.4 


I.Or 


0 . 8 - 


0 . 6 - 


CO 

*** 


0.4- 


0 . 2 - 


0.0 

c) 


A 


A 

4i 

d; 

£$j 

A/; 
A'* 


/. 


A 1 

A / 

a /; 

a / / 

A 

a 

A /' 

I ^ 


0 


I 

A 

t 

'a 

i: 

A 

£ 

A 


4 

4 

3 

A 

A|. 

a 1 ; 
A I. 




1 0 


I 

T 

T 

k 

b 

L 

t 

k 

A 

4 

4 

I 

4 

4 


A 


2. Or 


1 0 

u*/U 


i 






P/P 


\ \, 
\\ 

a: 


i 

( 

< 

k 

I 


A 
A 

l Si L 


12 3 4 


Figure 13. Comparisons of normalized velocity and pressure for cases R35D (blue long-dash line) and 
(green short-dash line) to experiment (open triangles). 


VII. Summary and Conclusions 

Several simulations were performed of a supersonic cavity flow exhibiting a reattaching free shear layer utilizing 
recent developments to the VULCAN CFD code. The simulations were designed in such a way as to: assess the 
applicability of LES/RAS to the chosen cavity flow, isolate the effect of inflow turbulence recycling on the flowfield 
statistics, determine any sensitivity of the solution to grid density, and evaluate the effect of shock sensor choice in 
the use of a hybridized flux reconstruction scheme. Two grids were used — a fine grid consisting of 35 million cells 
and a coarsened version of the same grid constructed of 8.8 million cells. A baseline run was performed on the fine 
grid using the Ducros sensor. Two cases were performed on the coarsened grid toggling the inflow turbulence 
recycling on and off in conjunction with the Ducros switch. The case performed on the coarse grid with recycling 
enabled was repeated using the Larsson switch to isolate the effect of shock sensor choice. 

Results that utilized the baseline fine grid compared quite well with experimental data outside of the cavity 
region. Through the cavity region, velocity profiles indicated good shear layer spreading rate agreement but a 
deviation in vertical location of the shear layer. This disagreement may be a result of the assumption of turbulent 
homogeneity in the spanwise z-direction or an insufficient spanwise width that could limit the size of spanwise 
structures. Even though measures were taken during the experiment to promote two-dimensionality, such as 
aerodynamic fences, previous studies have indicated that such approaches likely do not completely alleviate the 
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effects of three-dimensionality. For this case, it is easy to imagine three-dimensional motions arising from fluidic 
interactions with the comers of the finite-width cavity and with the fences. These same three-dimensional motions 
could potentially alter the dynamics of the shear layer. In future studies of this flow, three-dimensional effects 
should be considered to better capture these effects. 

Simulation results suggest that the solution is strongly grid-sensitive given the level of resolution utilized here. 
Hence, it cannot be confirmed that the finest grid employed is providing an acceptable level of grid convergence. 
However, simulations performed on the fine grid proved to resolve a large fraction of the turbulence, whereas the 
coarse grid simulations exhibited dominance of the subgrid modeled turbulence. Further, the coarse grid results 
suggest the inflow turbulence recycling plays a small role in the free shear layer development. Lastly, when used 
with the hybridized inviscid flux reconstmction method, the choice of shock sensor has little effect on the resulting 
flowfield. The sensors of Ducros and Larsson were considered for this work, and the results illustrate the relative 
insensitivity to the detailed functional form chosen for the shock sensor. 
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